//
// Author 
// Xianguo Lu
// lu@physi.uni-heidelberg.de 
//

#ifndef MATHBASE_H
#define MATHBASE_H

class MathBase
{
 public:
  typedef Double_t (*BFUNC)(const Double_t *, const Double_t *);
  typedef Double_t (*GFUNC)(const Double_t *);

  enum{
    NG = 2
    , NAG = 3
    , NAP0G = 4
    , NAP1G = 5

    , NAL = 3
    , NALG = 4
    , NALGE = 5

    , NMPVPUREION = 5
    , NMPVPURETR = 3
    , NMPVIONTR = 8
  };

  static Double_t G(const Double_t * xx, const Double_t * par);
  static Double_t AG(const Double_t * xx, const Double_t * par);
  static Double_t AP0G(const Double_t * xx, const Double_t * par);
  static Double_t AP1G(const Double_t * xx, const Double_t * par);

  static Double_t AL(const Double_t * xx, const Double_t * par);
  static Double_t ALG(const Double_t * xx, const Double_t * par);
  static Double_t DimlessALGE(const Double_t * xx, const Double_t * par);

  static Double_t MeanALGE(const Double_t * par, const Int_t nsample, TH1 *hh=0x0);

  static Double_t MPVPUREION(const Double_t * xx, const Double_t * par);
  static Double_t MPVIONTR(const Double_t * xx, const Double_t * par);
  static Double_t MPVPURETR(const Double_t * xx, const Double_t * par);

  static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);

  static Bool_t FindY(BFUNC func, const Double_t *par, const Double_t xmin, const Double_t xmax, const Double_t ytarget, Double_t &pos, Int_t &ndepth, const Double_t eps, const Double_t resolution);
  static Bool_t GetMPVFWHM(BFUNC func, const Double_t *par, const Double_t xmin, const Double_t xmax, Double_t *vout, const Double_t eps=1e-6, Double_t resolution=-999);

 private:
  static TString LocalName(){return "mathbasetmp";}

  //log(nor(1+eps)) = log(nor)+log(1+eps) = log(nor)+eps
  static const Double_t fEPS = 1e-6;
};

Double_t MathBase::G(const Double_t * xx, const Double_t * par)
{
  return TMath::Gaus(xx[0], par[0], par[1], kTRUE);
}

Double_t MathBase::AG(const Double_t * xx, const Double_t * par)
{
  return par[0]*TMath::Gaus(xx[0], par[1], par[2], kTRUE);
}

Double_t MathBase::AP0G(const Double_t * xx, const Double_t * par)
{
  return par[3]*( 1 + par[0]*TMath::Gaus(xx[0], par[1], par[2], kTRUE) );
}

Double_t MathBase::AP1G(const Double_t * xx, const Double_t * par)
{
  return par[3]*( 1 + par[4]*xx[0] + par[0]*TMath::Gaus(xx[0], par[1], par[2], kTRUE) );
}

//______________________________________________________________________________________
//_____________________________ Energy Loss Distribution _______________________________


Double_t MathBase::AL(const Double_t * xx, const Double_t * par)
{
  return par[0]*TMath::Landau(xx[0], par[1], par[2], kTRUE);
}

Double_t MathBase::ALG(const Double_t * xx, const Double_t * par)
{
  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*par[3]);
  const Double_t xmax = xx[0]+ng*par[3];

  TF1 ftmp(LocalName(),Form("[0]*TMath::Landau(x, [1], [2], 1)*TMath::Gaus(x, %.10f, [3], 1)", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, par, fEPS);
  
  return sum;
}

Double_t MathBase::DimlessALGE(const Double_t * xx, const Double_t * par)
{
  Double_t posp[NALGE];
  for(Int_t ii=0; ii<NALGE; ii++)
    posp[ii] = fabs(par[ii]);

  const Double_t gs = posp[3]*posp[1];
  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*gs);
  const Double_t xmax = xx[0]+ng*gs;

  TF1 ftmp(LocalName(),Form("[0]*TMath::Landau(x, [1], [2]*[1], 1)*TMath::Gaus(x, %.10f, [3]*[1], 1)*TMath::Exp(-[4]*x/[1])", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, posp, fEPS);
  
  return sum;
}

Double_t MathBase::MeanALGE(const Double_t * par, const Int_t nsample, TH1 *hh)
{
  const Double_t lanMean  = fabs(par[1]);
  const Double_t lanWidth = fabs(par[1]*par[2]);
  const Double_t gauWidth = fabs(par[1]*par[3]);
  const Double_t kapa = fabs(par[4]);

  TRandom3 rand(0);
  Int_t isample=0;
  Double_t sum = 0;

  while(isample<nsample){
    const Double_t de = rand.Landau(lanMean, lanWidth);

    const Double_t ew = rand.Rndm();
    if(ew > TMath::Exp(-kapa*de/lanMean))
      continue;

    const Double_t noise = rand.Gaus(0, gauWidth);
    const Double_t sig = de+noise;
    if(hh)
      hh->Fill(sig);

    sum += sig;
    isample++;
  }

  //TH1::GetMean is only accurate when under- and overflow are negligible
  return sum/nsample;
}

Double_t MathBase::MPVPUREION(const Double_t * xx, const Double_t * par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  //ALEPH parameterization: Double_t fPar[]={0.76176e-1, 10.632, 0.13279e-4, 1.8631, 1.9479};//ALICE TPC default
  const Double_t p0 = fabs(par[0]);
  const Double_t p1 = fabs(par[1]);
  const Double_t p2 = fabs(par[2]);
  const Double_t p3 = fabs(par[3]);
  const Double_t p4 = fabs(par[4]);

  const Double_t aa = TMath::Power(beta, p3);
  const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );

  //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);

  return (p1-aa-bb)*p0/aa;
}

Double_t MathBase::MPVPURETR(const Double_t * xx, const Double_t * par)
{
  const Double_t bg = xx[0];
  const Double_t gamma = sqrt(1+bg*bg);

  const Double_t p0 = fabs(par[0]);
  const Double_t p1 = fabs(par[1]);
  const Double_t p2 = fabs(par[2]);

  const Double_t zz = log(gamma);
  return p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
}

Double_t MathBase::MPVIONTR(const Double_t * xx, const Double_t * par)
{
  return MPVPUREION(xx, par) + MPVPURETR(xx, &(par[NMPVPUREION]));
}

//______________________________________________________________________________________

//Numerical Recipes in C, 1997, 3.1 Polynomial Interpolation and Extrapolation, p109
//Neville's algorithm
//error defined here depends on the choice of the route from ya[ns] to y_final
void MathBase::Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint)
{
  Double_t c[n];
  Double_t d[n];
  Double_t dift=-999;

  Int_t ns=0;
  Double_t dif = fabs(x-xa[ns]); 
  for (Int_t i=0; i<n;i++) {
    if ( (dift=fabs(x-xa[i])) < dif) { 
      ns=i; 
      dif=dift; 
    }
    c[i]=ya[i]; 
    d[i]=ya[i];

    if(kprint) printf("polint =====> x[%d]: %3.10e\ty[%d]: %3.10e\n", i, xa[i], i, ya[i]);
  }

  y = ya[ns];
  if(kprint) printf("======> y0: %3.10e at ns=%d\n", y, ns);

  Int_t did=ns-1;
  for (Int_t m=1;m<=n-1;m++){
    for (Int_t i=0; i<n-m; i++){
      Double_t ho=xa[i]-x; 
      Double_t hp=xa[i+m] -x;
      Double_t w = c[i+1]-d[i];
      Double_t den=ho-hp;
      if ( den == 0.0){
        printf("MathBase::Polint (den=ho-hp) == 0.0!!\n");
        exit(1);
      }

      c[i]=ho*w/den;
      d[i]=hp*w/den; 
      //if(kprint) printf("======> m %d, c[%d]: %5.6e\t d[%d]: %5.6e\n", m,i,c[i],i, d[i]);
    }
    if( m <= n -1-ns ){
      dy = c[ns];
      if(kprint) printf("------------> m %d, c[%d]: %5.6e\n", m, ns, dy);
    }
    else {
      dy = d[did];
      if(kprint) printf("------------> m %d, d[%d]: %5.6e\n", m, did, dy);
      did--;
    }
    y += dy;
    if(kprint) printf("======> y%d: %3.10e\n",m, y);
  }
}

Bool_t MathBase::FindY(BFUNC func, const Double_t *par, const Double_t xmin, const Double_t xmax, const Double_t ytarget, Double_t &pos, Int_t &ndepth, const Double_t eps, const Double_t resolution)
{
  //printf("MathBase::FindY %d: %.10e %.10e %.10e, %.10e, %e %.10e\n", ndepth, ytarget, xmin, xmax, pos, eps, resolution);
  if(ndepth>=50){
    printf("MathBase::FindY ndepth >= 50!! %d\n", ndepth);
    return kFALSE;
  }

  ndepth++;

  Bool_t kFindMax = kFALSE;
  Bool_t kFindMin = kFALSE;
  Bool_t kFindY   = kFALSE;
  if(ytarget>1e9)
    kFindMax = kTRUE;
  else if(ytarget<-1e9)
    kFindMin = kTRUE;
  else
    kFindY = kTRUE;

  const Double_t old = pos;

  Double_t x0 = xmin;
  Double_t dist =  1e10;
  Double_t ymax = -1e10;
  Double_t ymin =  1e10;
  Double_t b0 = -999, b1 = -999;
  while(x0<xmax){
    Double_t x1 = x0 + resolution;
    if(x1>xmax)
      x1 = xmax;

    const Double_t xevl = x0 + 0.5*(x1-x0);
    const Double_t yevl = func(&xevl, par);
    //printf("jjjjj %.10e %.10e %.10e %.20e\n", x0, x1, xevl, yevl);
    
    Bool_t kUpdateMax = kFALSE;
    Bool_t kUpdateMin = kFALSE;
    Bool_t kUpdateDist = kFALSE;
    if(yevl>ymax){
      kUpdateMax = kTRUE;
      ymax = yevl;
    }
    if(yevl<ymin){
      kUpdateMin = kTRUE;
      ymin = yevl;
    }
    if( fabs(yevl-ytarget)<dist ){
      kUpdateDist = kTRUE;
      dist = fabs(yevl-ytarget);
    }

    if( 
       (kFindMax && kUpdateMax)
       || (kFindMin && kUpdateMin)
       || (kFindY && kUpdateDist)
        ){
      pos = xevl;
      b0 = x0;
      b1 = x1;
    }
    x0 += resolution;
  }

  if(fabs(old-pos)>eps){
    return FindY(func, par, b0, b1, ytarget, pos, ndepth, eps, pow(resolution,2)/(xmax-xmin));
  }
  return kTRUE;
}

Bool_t MathBase::GetMPVFWHM(BFUNC func, const Double_t *par, const Double_t xmin, const Double_t xmax, Double_t *vout, const Double_t eps, Double_t resolution)
{
  if(resolution<0)
    resolution = (xmax-xmin)/100;

  Double_t ytarget = 1e10; //>1e9 for maximum
  Double_t pos = -999;
  Int_t ndepth = 0;
  if(!FindY(func, par, xmin, xmax, ytarget, pos, ndepth, eps, resolution))
    return kFALSE;

  const Double_t mpv = pos;

  //----------------- 

  const Double_t ymax = func(&mpv, par);
  ytarget = ymax/2;

  pos = -999;
  ndepth = 0;
  if(!FindY(func, par, xmin, mpv, ytarget, pos, ndepth, eps, resolution))
    return kFALSE;

  const Double_t leg0 = pos;

  pos = -999;
  ndepth = 0;
  if(!FindY(func, par, mpv, xmax, ytarget, pos, ndepth, eps, resolution))
    return kFALSE;

  const Double_t leg1 = pos;
  const Double_t fwhm = leg1-leg0;
  vout[0]=mpv;
  vout[1]=fwhm;
  vout[2]=leg0;
  vout[3]=leg1;

  return kTRUE;
}

//______________________________________________________________________________________
//______________________________________________________________________________________

#endif

/*
http://root.cern.ch/root/html528/TF1.html#TF1:Integral

Unless there is severe cancellation of positive and negative values of
   f(x) over the interval [A,B], the argument EPS may be considered as
   specifying a bound on the <I>relative</I> error of I in the case
   |I|&gt;1, and a bound on the absolute error in the case |I|&lt;1. More
   precisely
|G-I|/(Iabs+k)<EPS nearly always be true

TF1::Integral eps = 1e-1
ALG:    Real time 0:00:06.733, CP time 6.740
ALGE:   Real time 0:00:48.695, CP time 48.690

eps = 1e-12 (default in TF1)
ALG:    Real time 0:00:39.055, CP time 39.050
ALGE:   Real time 0:02:55.760, CP time 175.760

eps = 1e-6
ALG:    Real time 0:00:13.255, CP time 13.260
ALGE:   Real time 0:01:08.953, CP time 68.950

*/
